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S. Tanveer 1 

k Mathematics Department 

Virginia Polytechnic Institute & State University 

. Blacksburg, VA 24061 

Abstract 

An accurate analytic theory is presented for the velocity selection of a two dimensional 
needle crystal for arbitrary Peclet number for small values of the surface tension param- 
eter. The velocity selection is caused by the effect of transcendentally small terms which 
are determined by analytic continuation to the complex plane and analysis of nonlinear 
equations. 

The work supports the general conclusion of previous small Peclet number analytical 
results of other investigators, though there are some discrepancies in details. It also ad- 

I dresses questions raised by a recent investigator on the validity of selection theory owing 

to assumptions made on shape corrections at large distances from the tip. 
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1. Introduction 

The problem of the growth of a needle crystal in a pure undercooled liquid in the 
absence of any boundaries has received considerable attention in recent literature. In par- 
ticular, the growth of a steadily moving interface between solid and liquid has been studied 
using both analytical 1,2 ' 3 and numerical methods 4 ’ 5 . When surface tension is neglected, ex- 
act solutions with a parabolic crystal-melt interface were found earlier by Ivantsov 6 . These 
solutions form a degenerate set since for given undercooling and other experimental con- 
ditions, only the product of the tip radius of curvature and the steady dendrite velocity 
are determined in contradiction to experimental evidence 7,8 which suggests that each of 
these are separately determined for given undercooling far ahead of the interface. This 
degeneracy is not unexpected since in the absence of surface tension, there is not enough 
dimensional information to predict each of these physical quantities separately. 

When surface tension is taken into account, there is enough dimensional information to 
determine each of dendrite velocity and tip curvature in terms of undercooling. However, 
this need not imply that a solution exist in this case. Numerical evidence 4,8 appears to 
suggest that such solutions do not exist if we neglect the effect of crystalline anisotropy. 

Earlier, analytic study of phenomenological models 9,10 of solidification, suggested that 
solutions do not exist when anisotropy is neglected. The mathematical equations arising 
out of one of the phenomenological models 10 has been rigorously studied by Kruskal & 
Segur 11 . They prove that in the limit of zero surface tension, these model equations do 
not have any physically acceptable solutions when crystalline anisotropy is not taken into 
account even though the equations admit solutions when surface tension is exactly zero. 
This extraordinary situation arises due to the effect of terms beyond all orders in an 
asymptotic expansion for small surface tension. Kruskal & Segur extend earlier methods 12 
for linear equations to extract transcendentally small term in the asymptotic expansion of 
the solution to the third order nonlinear ordinary differential equation that they study and 
show that the leading order transcendental correction to a regular perturbation expansion 
fails to satisfy the condition on smoothness of the needle crystal right at the tip. However, 
when a term modelling crystalline anisotropy is included in the equations, a discrete set of 
solutions is found to exist. However, it is not clear to us that the simple model equations 
studied by Kruskal & Segur should faithfully reflect the properties of the actual needle 
crystal, even qualitatively. 

In the limit of small Peclet number, Pelce k Pomeau 13 reduce the original integro- 
differential equation called the Nash-Glicksman equation 14 to a simpler set of equations 
involving just one parameter. Subsequently analysis by Ben-Amar k Pomeau 1 of this 
equation and by Barbieri k Langer 2 of a simpler linearized form in the limit of small 
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values of a certain non-dimensional surface tension parameter support the conclusions 
of the numerical work at arbitrary Peclet number 4,5 for not too small surface tension 
that needle crystal solutions do not exist in 2 D or axi-symmetric 3-D case if crystalline 
anisotropy is neglected. (Numerical results become unreliable when surface tension is very 
small as the problem is nearly ill-posed in this limit.) On modelling the four fold crystalline 
anisotropy by a cosine term, the numerical work based on the Nash-Glicksman equation 
and analysis based on Pelce-Pomeau equations suggest that a discrete set of solutions 
exist for any nonzero crystalline anisotropy. Ben Amar-Pomeau’s analytical work formally 
extends the Kruskal-Segur 11 method for extracting transcendentally small terms to a 
non-linear integro-differential equation. This follows earlier work of Combescot et al 15 
who use the Kruskal-Segur method to the Saffman-Taylor finger problem, which again 
involves a similar non-linear integro-differential equation. The work of Barbieri et al 2 is 
based on an approximate linear equation and is based on Fredholm alternative condition 
on a non-homogeneous linear equation, where WKB approximate methods are used to 
find independent solutions to the homogeneous problem. This work follows the idea of 
Shraiman 16 , who employed a similar method for the Saffman-Taylor problem. Despite 
the apparent deficiency of such an approach in that the linear equations are approximate 
and that the WKB solutions are not quite correct in the neighborhood of turning points 
which must be encountered in evaluating the Fredholm condition using a steepest descent 
contour in the complex plane, the scalings in the dependence of physical quantities on each 
other turn out to be the same as the non-linear analysis of Ben-Amar & Pomeau, the only 
discrepancy being in the values of constants. 

However, for the axi-symmetric 3-D needle crystal, contradictory analytical evidence 
has recently been presented by Xu 17 . Rather than working with the Nash-Glicksman 
equation, he considers the original partial differential equations on both sides of the crystal- 
melt interface and obtains simplifications for small Peclet number using a slender body 
approximation. His analysis is not restricted to small surface tension. The basic approach 
used in his case is as follows: Given a slender axisymmetric 3-D crystal, he finds expression 
for the temperature in terms of an interfacial shape function. To the leading order, as Peclet 
number tends to zero, this expression is found to be a local function of the shape function, 
in contradiction with the Pelce-Pomeau simplification where the temperature at any point 
on the interface is expressed in terms on a global integral expression involving the shape 
function. Once the expression for the temperature on the needle boundary is found in terms 
of the shape function, Xu uses the Gibbs-Thompson condition of local equilibrium to reduce 
the problem to a 2nd order nonlinear ordinary differential equation which he then solves 
using phase plane analysis. He concludes that axi-symmetric 3-D needle crystals exist in 
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the absence of crystalline anisotropy and further, even for non-zero surface tension, each of 
dendrite velocity and tip radius is not a uniquely determined function of the undercooling. 
The degeneracy of the solutions is found to be the same as that for zero surface tension. Xu 
explains the discrepancy of his results with others by suggesting that the other researchers 
have implicitly assumed that the shape correction to Ivantsov parabola for nonzero surface 
tension tends to zero at infinity and thus have restricted the class of allowable shapes in 
their analysis and numerical computation. Thus the entire selection theory, atleast in axi- 
symmetric 3-D case, has been questioned. Despite some lively debate 18 , this controversy 
is yet to be settled conclusively. We find that Xu’s objections have some merit as far as 
the analytical evidence based on Pelce-Pomeau equation even for the simpler 2-D needle 
crystal problem, since in the derivation of Pelce-Pomeau equation from the original Nash- 
Glicksman equation, it appears to be necessary to assume that the shape modifications to 
Ivantsov solution tends to zero at infinity. Further, in the analytical work based on Pelce- 
Pomeau equation, the integrand is linearized based on the assumption that the correction 
to Ivantsov solution is small for sufficiently small surface tension. However, there does 
not seem to be any apriori reason to assume that for any small but fixed surface tension, 
the shape correction is small for the entire range of integration in the integral term and 
so linearization becomes questionable. As far as numerical evidence, Xu suggests that 
by truncating the infinite range of integration to a finite one and matching to the shape 
to a parabola at sufficiently large distances, one implicitly rules out shape corrections to 
the Ivantsov solutions that grow at large distances, though at a rate smaller that for a 
parabola. It is not clear to us if this argument has any merit or not. 

What is clear from all this is that one needs to resolve the discrepancy between the work 
of Xu and other researchers. Indeed, one can make a direct check on Xu’s leading order 
asymptotic expression for the temperature field on the 3-D axisymmetric needle boundary 
by a careful direct asymptotics of the integral term in the Nash-Glicksman equation for 
small Peclet number and checking if the expression is local or global. If Pelce-Pomeau ’s 
equation holds, then one needs to check the steps in Xu’s analysis leading up to the 
expression for the temperature field in terms of the shape function to find possible sources 
of error. This is currently under study. 

In the meantime, we thought it appropriate to reconsider the easier 2-D needle crystal 
problem, where Xu’s objections have some validity as well. Instead of considering Pelce- 
Pomeau simplification for small Peclet number, we thought it appropriate to consider the 
Nash-Glicksman equation for arbitrary Peclet number. When this paper was first written, 
we were unaware of any analytical work at arbitrary Peclet number, though the problem 
has been solved numerically. Since then we received a preprint of work at arbitrary Peclet 
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number by Barbieri & Langer 19 where they consider the needle crystal at arbitrary Peclet 
number in 2-D as well as 3-D using an approximate equation where the curvature term 
is linearized together with the integral term in the Nash-Glicksman equation. Aside from 
some quantitative errors in the value of constants that such a linearization would produce, 
their work does not address the objection of Xu as far as apriori assumption on the nature 
of shape correction at infinity. 

While our analysis is not mathematically rigorous either, we address some of the ques- 
tions regarding linearization of the integral term in the Nash-Glicksman equation with 
some care. Our final conclusions suggest that the selection theory, atleast in 2-D, is 
correct. The second order non-linear differential equation that arise in connection to cal- 
culating the leading order transcendentally small correction is found to be the about the 
same as that coming out of the Ben-Amar & Pomeau 1 analysis though their starting point 
was the simpler Pelce-Pomeau equation valid only for small Peclet number. We disagree 
with Ben-Amar &: Pomeau on several fine points in the analysis of this nonlinear equations 
particularly when the crystalline anisotropy is nonzero. In particular, the predicted quan- 
titative constant for the 1st branch of solution corresponding to dendrite with the largest 
velocity is found to be a little different from what we predict because we believe they use 
an analytic expression valid only for the higher branches of solution. 

2. Mathematical formulation 

In the frame of the steadily moving needle crystal, we fix the origin of the coordinate 
system (x, z ) = (0, 0) at the tip. The z axis be alligned in the direction of the crystal 
axis and the x axis perpendicular to it. A point on the needle boundary is descibed by 
the parametric representation (x(£), z(£)) , where £ is in the interval (— 00 , 00 ) and 

*(?) = e ( 1 ) 

z(0 = £ + x*(£) (2) 


This parametric representation is found to be rather suitable for avoiding nonuniformity 
in the linearization of Nash-Glicksman integral expression for the temperature as shall be 
seen shortly. The Nash-Glicksman equation for determination of the dendrite boundary 
can then be written as: 


A - <£o(l + a(l - cos 4 9)k = - f° d£'[l + ** { (0]« Pm '* U), ^o(P|r|) 

7T J- 00 


where 


r = (e-o 




1/2 


(3) 

(4) 
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In (3), 0 is the angle between the normal to the interface and the positive z direction 
and and k is the curvature of the interface (x(£), ^(£)) given by 


tan 0 = 


k — 


1 + 

1 + ~ £x fi gg 


(5) 


( 6 ) 


[1 + 2Xr£ + Xr\ + 

where the subscript with respect to £ denotes derivative with respect to £ . Note that 

8£ 2 (l + Xflg) 2 


1 — cos 40 = - r -» 

[l + 2 xr { + + £ 2 ] 

In equation (3), P (Peclet number) is defined as: 


(7) 


Ua 


2D 


( 8 ) 


where U is the velocity of the advancing interface, a the radius of curvature at the tip 
for the dendrite corresponding to the zero surface tension solution, D the temperature 
diffusion constant, do is the dimensional capillary parameter given by 


do — 


lcT n 

— TS 

Pa 


(9) 


where 7 is the surface tension, c the specific heat per unit volume; T m is the melting 
temperature, l the latent heat. In this paper, it will be assumed that for any given P , 
do is small enough so that do/P « 1 . This is not a severe restriction since the theory 
presented here is valid for small do . The parameter A in (1) is the non-dimensional 
undercooling defined as 

a = i(T m - rj ( 10 ) 


where T, » denotes the temperature at z = 00 far ahead of the finger. Note that each 
of x and z appearing in (3) are nondimensionalized by a . Also, in the definition of 
r in (4), the choice of a specific branch of the square root is made so that r > 0 for 
£ > £' and r < 0 otherwise for £ and £' on the real axis. Thus the absolute value 
|r| appearing in (3) is needed to be in accordance with the Nash-Glicksman derivation. 
The choice of a specific branch in (4) is made for the purposes of analytic continuation of 
(3) to the complex x plane as shall be seen later. 

When surface tension is neglected, i.e. d 0 = 0 , Ivantsov found exact solutions for 
a steadily growing dendrite with a parabolic interface shape with tip radius a (which 
is arbitrary). In our notation and non-dimensionalization, this corresponds to the exact 
solution 


Xr = 0 


( 11 ) 
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to (3). This is not immediately obvious on substitution of (11) into (3). However, Pelce 
& Pomeau have verified that (11) is indeed the solution to (3) provided the undercooling 
A is related to the Peclet number P by the relation 

A = n 1/2 P 1/2 e p erfc ( P l/i ) (12) 


Thus, when surface tension is neglected, it is clear from (8) and (12) that for given under- 
cooling and other experimental conditions, only the product of dendrite velocity U and 
the tip radius of curvature a are determined. However, experimental evidence 7 - 8 sug- 
gests that each of these two quantities are each separately determined as a function of the 
undercooling for other given experimental conditions. Thus, for an adequate theory, the 
degeneracy of these solutons needs to be removed. As pointed out earlier, any amount of 
surface tension introduces another parameter do into the problem and therefore there is 
then enough dimensional information for unique determination of each of U and a sep- 
arately. However, this does not guarantee that such a solution will exist and indeed our 
results suggest that in accordance with earlier numerical and analytical results (for re- 
stricted cases) , solutions exist only when the crystalline anisotropy parameter a ^ 0 . 

We now like to simplify the integral expression on the right hand side of (3) . We will 
assume that for small d 0 and fixed £ , xr is small. However, as shown in the appendix, 
the boundary condition that the non-dimensional temperature on the interface approach 
A , a constant, as as £ — ► ±oo can allow for the interface shape correction function 
xr to grow with £ at a rate like £ 1- ' , where s > 0 . Thus, xr , need not uniformly 
be small. However, it is reasonable to assume that xr^ is small uniformly for all £ , and 
thus from expression (4) for r on application of the mean- value theorem on the quantity 
5 it is clear that the deviation from r from ro is small for small do for all 
£ and £' , where 


ro = («-£') 


i + j(£ + £')* 



(13) 


Thus it is legitimate to linearize the right hand side of (3) for any given P for sufficiently 
small do . If we subtract off the Ivantsov solution, we find that to linear order in x R on 
the right hand side: 


-d 0 [l + a(l - cos A6))k = - f° dt'xR'[Z')e p W'>- z WK 0 {P\ro\) 

7T J-oo 

+ —T d(’ (£ - OM£') - **(£)] (14) 

7T J-oo |ro| 

Note that in (14), the subscript with respect to £ denotes derivative with respect to £. 
In (14), we used the identity that derivative of K 0 is -Ki . It is convenient to get rid 
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of the x R £ term by integrating by parts and we find after careful consideration of the 
singular nature of the integrand that (14) is equivalent to 


— d 0 [l + «(1 - cos 40)]/c = x*(£) — [°° d^'e p ^~ z ^ (£' - £) 

7T J-oo |ro( 


- ^ f«-o(F|r 0 |) + gl( |f o | rol) [»(e) - *«’)] 


(15) 


Now for the Ivantsov solution, the nondimensional temperature within the crystal is a 
constant, A and so 

A = - f°° d? e p W')- 2 V Ko(P\f 0 \) (16) 

7T J —oo 

where (x, z ) is now inside the crystal and 

,1/2 


r 0 = 


(*-« ! + (*+j« ,! ) 


The partial derivative of the above expression with respect to x must be zero, since the 
temperature within the crystal is uniform for the Ivantsov solution, when the curvature 
effects are neglected. On the other hand if we take the derivative of the right hand side of 
(16) with respect to x and approach the interface from the inside of the crystal we find 
that 

* |r„| f, “ £ + (1 + z\) “ ° 

Thus (15) can be further simplified as 


(17) 


-d 0 [l + o(l — cos 40 )\k = ““ 2 rZ|i(£) — — / 

[1 + Z c ] 7T J-oo 




(i + 

xoirkoi) + m ^te) - *(e')] 


d(' 


(18) 


3. Regular perturbation expansion and analytical continuation to the upper 
half £ plane 

If we now carry out a regular perturbation expansion of xr in powers of do : 


xji(£) = doZiU) + d 2 0 x 2 {t) + ... 
we find that Xi satisfies the linear singular integral equation 

-[1 + a(l - cos 40 o )]ko = £ 

Kx(P\ro\) 


(19) 






tfo^M) + 


ko| 


we) - 2(01 


d? 


( 20 ) 
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where 8 0 and k q are equal to the expressions (5) and (6) for 8 and k with the 
substitution xr — 0 . 

We numerically calculated a smooth solution to (20) by discretization and satisfying 
the equation at a discrete set of points. In addition to (20), we imposed the condition 
ii(0) = 0 so that the tip of the dendrite coincides with ( x,z ) = (0,0) . The resulting 
linear system was solved without any difficulties and consistency of the solution checked 
by doubling the each of the number of discretization points and the size of the truncated 
domain. The solution, as expected, was found to be an odd function of £ implying a 
smooth symmetric dendrite at least to order do . In particular, this implies that the tip 
of the needle crystal is smooth. It is conceivable that the same is true to every order in 
the expansion (19) though we have not calculated higher order solutions. We assume that 
this is indeed the case. 

At this point, it is appropriate to point out that if instead of the parametric represen- 
tation (x(£), z(£)) , for the free boundary, we had used (x, z(x)) representation and 
decomposed 

*(*) = x2 + M*) ( 21 ) 

and carried out a linearization of the integral term in the Nash Glicksman equation, we 
would arrive precisely at (18) with —Z(Xr(Q replaced by zr(x) , with £ and £' replaced 
by x and x' and z(£) replaced by — |x 2 . However, justification of the linearization of 
the integral term of the Nash-Glicksman equation for such a representation appears to be 
difficult if such a representation were used. 

Now, we proceed to calculate the leading order transcendentally small correction to 
(19). Following the ideas of Kruskal & Segur, we do so by analytically continuing (18) 
to the upper half x -plane to find sources of nonuniformity of the expansion (19). These 
sources of nonuniformity in the complex x -plane contribute transcendentally small terms 
in the asymptotic expansion of Xr and it is our intention to calculate the leading order 
transcendentally small term in order to find any constraint on the parameter d 0 arising 
from the requirement that the tip of the parabola be smooth. It is convenient to define 

2je(£) = £ £*(£) ( 22 ) 

Note that Zr(£) is not defined as z#(x(£)) ; however when £ = 0(1) , to the leading 
order in d 0 , the two are the same. Note that there can be deviation of z R (£) from 
Zr(x(£)) which is not uniformly small for all £ even for small do . 

If we restrict our attention to symmetric needle crystals for which zr(£) = zr(— £) and 
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substitute (22) into (18), then 


p 2 r° i, 


+ 


1 f ifi(F|r 0 |) 

2 1 kol r, | IS 


where 


-«' 2 1^(£') [JT 0 (P|ro|) + K 0 {P ri ) 

-e\ 

<*£' 

(23) 

»i i/2 





(24) 


We note that with the choice of branch in the above squareroot ri > 0 , for £ and 
£' on the negative real axis. Further 


r = r 


and for £ < £' , we have to choose 


r = r e 


(24) 


(25) 


in order to analytically continue to the upper half complex £ -plane. It is well known that 

K 0 {r e~ ir ) = K 0 {r) + i tt I 0 (r) (26) 

Thus (23) can be seen as the limit of £ approaching the negative real axis from the upper 
half complex £ plane of the following equation 

— d 0 [l + <*(1 — cos 40)]/c(0) = /(£)^(£) + f° d£'G(£, £>*(£') + f° d? J(t, £')^(£') 

(27) 

where 


G(f. (') = — fit„(Fr„) + MPr I0 ) + i ( _ 

it [ 2 ( ro r i J 


e 2 


and 


no = 


i — if 


(28) 

(29) 

(30) 


From (20), in an analogous procedure, it is found that the analytical continuation of 
the leading order regular perturbation solution zi(£) = £xi(£) in the upper half complex 
£ plane satisfies 

-|l + a(l - cos 4» 0 )|<O = *(f)*i(f) + f ((') + /V J(0 «')*.(«') (31) 

J -OO J t 
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It is easy to see from (31) that Zi(£) is singular at £ = t in the upper half complex 
£ -plane. From symmetry of the equation, it is easy to see that z\ is also singular at 
the lower half complex £ plane at £ = — * . Thus, we need to find local equations in 
the neighborhoods of these points such that as the real axis is approached, the solution 
matches with the regular perturbation expansion (19). The terms that will not match 
must be transcenden tally small in the physical domain. 

To find the form of the leading order transcendentally small term, we subtract d 0 times 
(31) from (27) assuming that zr is a small deviation from d 0 z x to find that the resulting 
homogeneous part of the equation for the small deviation zg is 




1 + 


8a£ 2 


(i + ey 


+ 4 


H 


3 £ 


56a: £ 3 16ft£ 

/, i 


+ 


= Z H 


(l + £ 2 ) (l + £ 2 ) 3 (l + £ 2 ) 2 J d 0 (l-it) 

3 56ft£ 2 16ft 


where 


and 


(i + £ 2 ) (i + £ 2 ) 3 (i + f 2 ) 2 

ll + £ 2 1 3 / 2 ro 

= - <?(£,£») **(0 
CLq J - oo 


+ h + h 


(1 + 

(32) 

(33) 


h = 


r m' j(W) z„a') 

J £ 


(34) 


_ (i + a s / * 

d o 

4. Transcendentally small terms for a = 0 

For ft = 0 , i.e. no crystalline anisotropy, the leading order asymptotic solution for 
small d 0 to the linear integro-differential equation (32) in the upper complex £ plane away 
from the immediate neighborhood of the turning point £ = t must be linear combinations 
of and 02 defined as: 

9ll = + (35) 


Note that the above is just the two independent WKB solutions to (32) with the right 
hand side of (32) neglected. On substitution of (35) back into 7 4 and h it is clear that 
these contributions are of smaller order in do compared to other terms on the left hand 
side of (32). We note that on the imaginary £ axis in the interval (0,t) 0 2 is real and 

transcendentally small, while 0 i is transcendentally large. This is also true for a certain 
region in the complex £ plane in the neighborhood of the imaginary axis in the interval 
( 0 ,t) (sector I as sketched in Fig. 1 ). However, there is another region in the complex 
£ plane with Re £ < 0 , where 0 i is transcendentally small and 0 2 transcendentally 
large (sector II in Fig. 1 ). The boundary between these two sectors is called a Stokes line 
and is determined by the condition 

Re {* J* (1 + £ ,2 ) 3/4 (1 - t£T 1/2 > = 0 ( 36 ) 
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The sketch in Fig. 1 is justified from the following consideration: First we note that 

*■ / V ( i + «' 2 ) s/4 ( 

is purely real and positive. Again, consider the real part of 

Rc{i (i -i- (i - .i') -17 * } < o 
Jo 

for £ , on the negative real axis. By considering the argument of (1 — i£') -1 / 2 } one easily 
establishes that the left hand side of the above equation it is a monotically decreasing 
function of |£| on the negative real axis. Using this, it is easy to show that only one of the 
Stokes line emanating from the turning point £ = t intersects the negative real £ axis 
as shown in Fig. 1. 

Including the leading order transcendentally small correction, in sector I (that includes 
the imaginary £ axis between 0 and i ), 


zr ~ d 0 zi + HOCP + Ci g 2 (37) 

where HOCP stands for higher order corrections with power dependence in do . From 
now on, we do not bother to write HOCP though such terms are present in the expression 
for zr and dominate the leading order transcendental correction in d 0 which will be 
explicitly written down as they determine the velocity selection. 

As £ = 0 is approached, the importance of the transcendental term arising due 
to the effect of singularity at £ = —i becomes as important. At exactly £ = 0 , the 
contribution of singularities at £ = ±* are of the same order in d 0 . Since z R must be 
real on the entire real £ axis, it follows that on the real £ axis in some neighborhood of 
£ = 0 , the contribution from £ = —i to the leading order must be C{ g 2 ( * denotes 
complex conjugate) so that on the real £ axis 


Zr ~ doZ\ + C 2 </2 + C 2 g 2 


(38) 


is real. It is easily seen that the slope at the tip as we approach it from the negative £ side 
is 


dz R 

di 


2Re [Ci ri(0)] 


On using (35) 


_ 2 Im C l P 1 ' 2 V /! 


(39) 


Thus a smooth tip implies 


Im Ci = 0 


(40) 
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In sector II of Fig. 1, which borders on the real £ axis for sufficiently negative £ , 
including the leading order transcendental small correction, we must have 

Zr ~ do Zi 4- C 2 <7i (41) 


and in this sector, on the real axis we must have 


zr ~ do Z\ + C 2 <71 + C\ g\ 


(42) 


To find Ci and determine if the smooth tip condition (40) can be satisfied, we must 
consider the immediate neighborhood of £ = » in the upper half plane, where each 
of the expressions (37) and (41) are invalid both because of the linearization used in 
obtaining (32) and the fact that £ = t is a turning point. We introduce local dependent 
and independent variables variables F and f defined as 

{ = .’(I - 4 ,T P-V 2-*/' S ) (43) 

zr = -d' 1 ’ P" </ ’ 2~ !/i F (44) 

Then it is found that (35) is to the leading order in d 0 reduced to 

F" - (f - F') S/2 F - 1 (45) 


In obtaining (45) from (35), the contribution from I 4 is of order d 0 since I 4 involves 
the integral of z R on the real axis where zr = O(d 0 ) . As far as 7 5 , one needs to 
be more careful, since in the range of integration including the immediate neighborhood 
of £ = * where the scaling (44) holds. However, on carefully analysis, it is found that 
I5 does not contribute anything to the leading order as well. It is easily seen that the 
asymptotic behavior for large f that matches with zr = do z\ is 


F ~ 


1 

^3/2 


(46) 


To find transcendentally small correction to this, we linearize (45) about (46) and find that 
the homogeneous part of the linear equation is 

Fh" - y ( F h ’ ~ ? ,2F » = 0 ( 47 > 


The transcendental correction to (46) to the leading order for large f must be linear 
combinations of the WKB solutions to (47) given by 


f 3/8 e ±±fV4 


(48) 


13 


For Arg f in the interval (— 2 tt/ 7, 0] , inclusion of the leading order transcendental 
correction gives 

(49) 


C 3/! 


L- + A lS >l‘ 


This matches with (37) in sector I provided 


_ _ 2~ 6 / 7 p - 13/28 ^ 13/28 


(50) 


For large f with Arg g in (— |rr, |jr) , 

* ~ -^ + x 25 ’/ 8e » t - 
and this matches with (41) in sector II of Fig. 1 provided 

= _ 2 ~ 6/7 p - 13/28 ^ 13 / 28 

4 ; ° 


(51) 


(52) 


We note that we are interested in a solution to (45) which for large f with Arg f in 
(— 6 it/1, 0]) has the asymptotic behavior given by (46). Our numerical calculation of 
appropriate solution to (45) involved solving (45) on two rays emanating for f = 0 going 
to a large distance from the origin with Arg $ = -tt/2 and Arg $ = 0 with asymptotic 
boundary condition (46) at the other end points of these straightline contours. For a given 
trial value of F’(O) the two point boundary problem on each ray was solved by standard 
second order discretization of (45) and using Newton iteration. Once these solutions were 
obtained numerically, one sided second order differencing gave us the estimated value of 
F'( 0) on each of the rays. In an outer Newton iterative procedure, the trial value of 
F( 0) was found so that the computed F'(0) along the two rays agree. From monitoring 
the size of the Jacobian, it was clear that the problem was not underdetermined and we 
checked that indeed a unique solution to (45) satisfying given decay conditions exist. Once 
the solution converge, the imaginary part of solution F along the ray coinciding with the 
positive real f axis at large distances was found to proportional to $ 3 / 8 e~ $ {7/ * with the 
proportionality constant equalling —0.875 . From (49), it follows that Im A x = —0.875. 
From (48) and (50), the tip slope 


^2(0) = — 2 1 ' 7 Im (A x ) V p ' n S° 

d£ 

which on numerical evaluation is 


= 2 ,/7 P'/ 2 %' ,,a 0.875 e" 06156677 d °''' 


(53) 


which is clearly non-zero. Thus no needle crystals exist since a jump in the slope at the tip 
implies infinite curvature which means that (l) could not possibly be satisfied at the tip. 
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The formal solution that we have constructed is an asymptotic solution of (l) for £ real 
in the interval (— oo,0) where we relax the requirement of a smooth tip. The same result 
with almost the same numerical values was obtained by Ben Amar &: Pomeau for small 
Peclet number. Here, we see that (53) holds even for arbitrary Peclet number. 

However, we differ with Ben Amar &; Pomeau’s analysis on a certain point which does 
not change the result the result (53) but is important as far as checking consistency of 
solution. They claim that the solution F to (45) is singular when f — ► 0 and find 
the need of an inner neighborhood with a different scaling. Their argument is based on a 
possible behavior of (45) near the origin. However, not every solution to (45) need have a 
singular behavior at the origin and indeed from numerical integration of (45) (with careful 
choice of consistent branch cut) , we find that the solution to (45) that satisfies the decay 
conditions at oo for Arg f in [0, 67r/7) remains finite at f = 0 . Indeed, if F tends 
to oo as f — * 0 , the linearization of the integral term in Nash-Glicksman equation or 
even the Pelce-Pomeau equation for complex £ in the neighborhood of £ = t would 
then be questionable. 

However, this discrepancy with Ben-Amar Pomeau’s analysis has no bearing on the 
final result (53) which are in agreement. 

5. Transcendentally small correction for nonzero anisotropy 

The WKB solutions to (32) for small do are now given by g\ and g 2 , where 

01,2 = (1 + e 2 ) 7/4 £' 1/4 e~ 1/2f < d *' m ' )+ A ] e ±( d o 1/a />' L ' /3 M (54) 


where 


L = 


1 - t£ 


- (i + e 2 ) 3/2 i + 


8a£ 2 

(i + e 2 ) 2 J 


i-i 


( V 

56a£ 3 \ 

fl + 8 “£ 2 1 

l l + ( 2 

(i + e)v 

1 (i + e 2 n 


(55) 


(56) 


Note that each of L and Q are singular at £ = £ 0 on the imaginary £ -axis between 0 
and i , where 

(57) 


fo = •[(l+2a) ,/J -(2a) 1/J ] 


The WKB solutions are invalid in a small neighborhood of £ = t and £ = £o • The 
form of the local equations depend on the size of a . In the next two sections we consider 
two cases: a P 4 / 7 d^ 7 = 0(1) and a P 4 / 7 d^ 7 » 1 : 

6. Transcendental correction for a P*P do* /7 = 0(1) 

In this case £ 0 is within a d^ 7 neighborhood of i . The the WKB solutions (54) 
holds beyond a d^ 7 neighborhood of £ = * as in the previous section. To the order 
of approximation to which (54) is valid, we can replace g\ and g 2 given by (54) by the 
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simpler WKB solutions gi and g 2 as in (35). This is because in (55), a is small and 
the terms involving a are only important near £ = t , where the WKB solutions are 
invalid any way. Near £ = t , we introduce the same change of variables (43) and (44) 
to find that the leading order equation is now 


F" - 


(£ - f ) T/ * 

l(r - n ! - P\ 


F = 1 


(58) 


where 


/ 3 = 2 9 / 7 a P 4 / 7 do 4/7 (59) 

As before, the asymptotic behavior of (58) that matches with d Q z\ when d%j 7 « |1 + 
*£| << 1 is 

(60) 




We linearize (58) about this solution and obtain the transcendental correction to (60) from 
the WKB solutions of the form (48). Once again, as in the previous section, (49) is valid 
for large g with Arg g in (— 27r/7, 0] , and this matches with 


Zr ~ do Z\ + Ci <72 (61) 

in sector I (Fig. 1) provided (50) holds. Similarly, for large g with Arg g in the 
interval (— 67r/7, — 27 t/ 7) , (51) holds and this matches with (41) provided (52) holds. 
Thus a unique solution to (58) is found by requiring that the solution goes to zero for large 
g with Arg g in the interval (— 6 tt/ 7, 0]. However, for such a solution for arbitrary /? , 
we generally have Im A x ^ 0 implying Im C\ ^ 0 . This implies a non smooth tip 
in the general case. However, on varying /? , we obtain a set of values of (3 and hence 
d 0 for given a for which the smooth tip condition is satisfied. The smallest /? value were 
found numerically to be 1.4926. The details of the numerical method is given in section 7. 
Note that the scaling of do with a follows from the definition of (3 and is consistent 
with earlier numerical 4,5 and analytical work 1,2 . The results for the case of large (3 with 
a << 1 is a special case of the case considered in the following section though it can be 
treated by a direct analysis of (58). 

7. Transcendental correction for a P 4 / 7 d^ ^ 7 » 1 

Note that in this case, we could either have a = 0(1) or a « 1 provided /? as 
defined by (59) is very much larger than unity. At the outset, we will be assuming that 
a is order unity. Later, scrutiny of the assumptions show that the final result is valid even 
for small a provided (3 is large. 
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In this case, the WKB solutions (54) do not simplify to (35). The Stokes lines are 
determined by the condition 


Re i J* d? 


= 0 


where L is given by (55). The Stokes lines in this case are shown Fig. 2 and the asymptotic 
growth shown in sectors I and II are now relevant since they extend all the way to the 
negative real £ -axis. 

In this case, we introduce the independent and dependent variables in the neighborhood 
of £ = i given by 

« = «•(!- 4 /u 2 1 ' 11 P~ lln a */u C ) (62) 

zr = -4 /u 2 ! ' u P- , I" a*'" F (63) 

Then the leading order equation for £ of 0(1) is 


F" + (f — F') Vi F = 1 


(64) 


For large f , the asymptotic behavior that matches with zr ~ doZi when 4 /T « 
|1 + *£| << 1 is 

F ~ (“) 


To find transcendentally small corrections to this behavior, we linearize (65) about this 
leading order behavior and find WKB solutions to the homogeneous 2nd order linear 
ODE. Including this transcendental correction, we find that for large £ with Arg $ in 
(-47T/11, 0) 

F ~ + A, (66) 


and this matches with 


in sector I of Fig. 2 provided 


f 7/2 

Zr ~ d 0 Zi + Ci 02 


— = pi/M 2 ~ 39 / 22 a " 1 / 22 do /44 

-Al 


(67) 


( 68 ) 


For large f , for Arg $ in the interval (— 87r/ll, — 47T / 11) the leading order behavior of 
F is given by 

F ~ -777 + A 2 <r 7/8 e‘^ fU/4 (69) 


and this matches with 


f’/! 

Z R ~ ^qZ\ + C; (ft 


(70) 
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in sector II (Fig. 2) provided 


£* = - e ’ v / 4 P 1 / 22 2" 39 / 22 a " 1 / 22 4 /44 (71) 

■^2 

Thus a unique solution to (64) is calculated by requiring that the asymptotic behavior of 
the solution F for large <; be given by (65) with only transcendentally small correction 
for Arg f in the entire interval (— 87r/ll,0) . It is clear that A\ determined as such can 
only be a pure number. We do not determine this pure constant A\ in this paper. 

From the arguments similar to that of section 4 leading up to (40), it clear that the 
appropriate condition for a smooth tip is that (67) be real on the imaginary £ axis near 
£ = 0 . Thus it is necessary that 


Arg 


\ Cx e - ido ~ 1/3 //“ L ' /2 (€') 


— n 7 r 


(72) 


where n is some positive integer. Note that the choice of a negative sign on the right 
hand side of (72) follows from the sign of 


f '° dC L'I'(C) 


which is negative because with the choice of branch, Arg L 1 / 2 varies continuously from 
0 to 7r/2 as Arg (£ — £ 0 ) varies from 0 to — ir . Thus, from (68), we find that the 
condition of smooth tip implies that 


where 


where 


P 1/2 do 1/2 G(a) = n 7r + ^ + Arg A\ 

4 


G(a) = f dy 
Jyo 


(i - »*) s/< 

d +*) ,/4 


8 ay 2 

.(i - y 2 ) 2 



Vo = £oA = yO- + 2a ) - >/2a 


(73) 

(74) 

(75) 


Equation (73) is the selection rule. 

Equation (73) is also valid small a only if /? as defined in (59) is very much larger 
than unity. In the cause for small a , it is easily seen that 


G(«) ~ 2 9 ' 8 o'/* £ dq 
which on numerical evaluation gives 

G(a) - 1.80205a 7 / 8 

I 

I 
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(76) 


(77) 


I 



In terms of /? , (73) reduces to 

7T 8/7 „ 

/? = 1.2437(n7r + — + Arg AJ (78) 

4 

Ben-Amar & Pomeau also arrive at the result (73), but they implicitly assume that 
(73) is valid for any a and for any value of integer n . We claim that (73) can only 
strictly hold for large values of n because if n were of order unity, G(a) in (73) will 
have to be small and of order d\j 2 P~ l ! 2 (which has to be small for the theory to be valid). 
From (76), this would imply that /? = 0(1) and then the result (73) is not strictly 
valid. In this case, one has to use the results of section (6) provided a is small. If a is 
not small, we cannot use any of the results of this paper or the previous ones 1,2 to make 
a proper prediction of do for the 1st few branches of solution i.e. n = 0(1) , since 
the corresponding do P~ l are not small and therefore beyond the validity of the theory. 
However, despite the fact that (73) is strictly invalid for n not large, it appears from 
comparison with direct numerical calculations 1 that the formula is surprisingly accurate 
even for relative small n and a over the range of experimental conditions. 

Notice that the asymptotic form of solution (67) is also invalid near f = fo as pointed 
earlier by Ben Amar & Pomeau. However, this point has no bearing on the result (73). If 
we are interested to find the behavior of the solution in this neighborhood, we introduce 
local variables 



t = £o - itd 0 x 

(79) 

and 

z R = -t 2 d 2 0 G{x) 

(80) 

where 

-i6*ae 0 (i-eg)(i-te 0 ) 
.p(i + (if* 

(81) 

Then to the leading order in 

d 0 , the equation for x = 0(1) is 



(G" - 1) (x " G') = G 

(82) 


For large x > G ~ — X and linearizing (82) about this and finding WKB solutions to 
the associated homogeneous equation, we arrive at the following expression for G for 
large x that (with approriate choice of constant B ) matches with (67) as £ — *■ £o f° r 
Arg x in (-», 0] : 

G x + B x _1/4 e" 2 * 1 ' 2 (83) 

This does not affect the selection rule (73). 

8. Numerical determination of /? of order unity 
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Here in this section, we describe the numerical method used to determine /3 so that 
solution to equation (58) satisfies the asymptotic condition (60) for large f for Arg f in 
(— 67r/7, 0] and that the solution be real on the positive real f axis for f sufficiently 
large. From Schwarz reflection principle, it follows that we are interested in a solution that 
satisfies the asymptotic decay condition (60) for large f when Arg $ in ( — 67 t/7, 67t/ 7) . 
Since only one such solution could be found, it follows that solution satisfying decay con- 
dition (60) for Arg f in (— 67r/7, 6 tt/ 7) must automatically satisfy condition that F be 
real for sufficiently large f on the real axis and indeed that was checked numerically. 

The method employed is similar to the one employed earlier 20 in the context of the 
Saffman-Taylor finger problem. We choose a point xo on the positive real axis that is 
sufficiently large so that the resulting solution is real at f = z 0 . This was done by trial 
and error. However, we do not choose xq unnecessarily large because such a choice will 
cause numerical inaccuracy. 

We go through the procedure given in the next two paragraphs to calculate the residual 
corresponding to a given value of /? : 

We take N points lined up parallel to the Im f axis of the form ft = Xq — iL\ + ikh , 
where k is an integer ranging from 0 to N + 1 , Li is a large positive number far larger 
than x 0 , and h is the distance between adjacent ft points. N is chosen to be an odd 
integer and h chosen so that (N+l)h = 2Li . The asymptotic condition (60) is employed 
at the end points ft and ftv + i and the (58) discretized and satisfied at f = ft for 
k ranging from 1 to N using standard second order finite differencing. This discretized 
two point boundary value problem is then solved using Newton iteration choosing an initial 
guess F — 0 and convergence was obtained without any problems. Once convergence is 
attained, we store the value of F and its estimated derivative obtained by second order 
central differencing at f = z 0 • 

The same procedure as in the last paragraph was used for a set of points on the real 
axis, ft = z 0 + jhi , where j now ranges from 0 to Ni + 1 , with (Ni + l)hi = L 2 , 
where is a large positive number and N\ is a large positive integer so that hi is small. 
The equation (58) is discretized and satisfied for j = 1, ...Ni and the decay condition 
(60) used at end point corresponding to j = N\ + 1 . At j = 0 end point, we use 
value of F as obtained in the last paragraph. Once a converged solution is obtained 
on this contour, we estimate the derivative of F at f = xo by a one sided second 
order differencing. The real part of the difference of estimated derivative here and in the 
procedure of the last paragraph is the residual. The imaginary part is automatically zero 
to within machine precision, as it must be from the symmetry of the equation. 

Once the residual is calculated for given (3 , in a Newton iterative procedure, we drive 
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the residual to zero. The smallest value of /3 so found was 1.4926 and this corresponds to 
the dendrite moving with the largest velocity. We took N and Ni to be 2049, Li and 
L 2 to be each 10 and x 0 = 2.0 and the results were unaffected by doubling each of N , 
Nx or by changing Li , L 2 and x 0 . We do not carry the calculation for other branches 
because experience has shown that (73) becomes quite accurate even for moderate values 
of n though the expression should only be asymptotically valid for large n . 

9. Discussion and Conclusion 

We present here an analytic theory for the determination of velocity for two dimensional 
dendrite at arbitrary Peclet number in the limit of small values of the surface tension 
parameter provided the ratio of surface tension and Peclet number is also small. We point 
out some discrepancies with earlier analytical work carried out in the limit of small Peclet 
number. The method is both qualitatively and quantitatively accurate and is an attempt to 
answer some serious objections raised by an earlier investigator on the validity of selection 
theory. 
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Table 1 


a G(a) 

0.00001 0.0000756 
0.0001 0.000138 

0.001 0.00404 

0.01 0.0266 

0.02 0.0451 

0.03 0.0605 

0.04 0.0738 

0.05 0.0855 

0.06 0.0961 

0.07 0.106 

0.08 0.115 

0.09 0.123 

0.10 0.130 

0.12 0.144 

0.14 0.155 

0.16 0.166 

0.18 0.175 

0.20 0.183 
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Figure 1: Solid Lines: Stokes Lines determined by Re{i // d£'(l + £ f2 ) 3 / 4 (l — »£') 1 / 2 } 

= 0 in the complex £-plane for Re£ < 0. 



ed by Re[iff L 1/2 (^)d^'] 


Define 


Appendix I 


«[fl = 7Jl dli 'l 1 + W ]e '' rW ‘~ e) Jr ° (p ' r ') 


dxj 

dt 


where 


<• = [(£ - (f + 2 (f - 0 MO - **(«') + (**(£) - **(£'))’ + \ (0 - £’V 

The purpose of this appendix is to show that the boundary condition 

lim u(£) = A 

£ — ► ±oo 


(A 1) 

• 1/2 
A2) 

(A3) 


, can be satisfied for any shape correction function £Er(£) satisfying the following bounds: 

MOI < c, [i + lei 1 -] ( m ) 


and 

dx R ^ C 2 

dt 1 i + ier 

for some constants C\, C 2 and s with C\ > 0 , 1 > C 2 > 0 and | > s > 0 . 
Note that the A appearing on the right hand side of (A3) is related to P through (12). 
Thus the shape correction from the Ivantsov parabola can actually grow at oo . The 
upper bound on C 2 is not too restrictive since the function is expected to be small 
for d 0 reasonably small, though the proof does not assume anything directly about the 
size of d 0 . 

We will carry out the proof only for the symmetric dendrite, .i.e. when xr(— £) = 
—Xr(£) though it is true in general. For a symmetric dendrite (Al) reduces to 

«[£] = - 1° d(' [1 + e-W’-?) [«,(P| r |) + if 0 (p ri )] (.46) 

7 T J — oo CiQ 

where 

r\ = [(£ + £') 2 + 2 (^ + CO (*Jt(0 + + (**(£) + *Jt(0)* + T (£ 2 “ C' 2 ) 

(A7) 

Choose any Si satisfying the condition 



1 > Si > 1 — s 


(AS) 
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and choose 

m = (-er m 

We assume — £ is large enough so that — £ + n < 0. We now decompose the function 


u 


where 


«ll fl = 7 <*€' [1 + ^f\ |Ko(P|r|) + Jf„(Pr,)] (All) 


U = tti + «2 + Us 


dx R . 


(-410) 


dt 


« 2 [ t] s - f* +tl d£' [1 + e-iW 3 -^ [K 0 {P\r\) + K 0 (P ri )} (A12) 

M - 7 Ji di ' I 1 + ^1 [Sb(P|rD + K-o(Pri)] 


dt 

dXR y _ 

dt 

Now, from the properties of Kq it is clear that there exists constant B\ such that 


(A13) 


0 < P 1 /* |r| 1/J K 0 (P|r|) c p ^ < fl, 
0 < P 1/! r, 1 ' 2 Ko(Pr t ) e r, ‘ < B l 


(A14) 

(A15) 


Thus 


p 1 / 2 ^ p€-m e - P (( 

|«i| < (1 - C 2 ) B X 2V2 / d£' *-= 

7T J — oo * 


.-p(C’-e’) 


\[F - « ! 


< (1 - C,) P‘/ ! B, 2 ^ /” dg 


e- r i 

y/i 


(A16) 


It is clear that, in the limit of n — ► — oo , the right hand side of (A16) goes to zero. 
Now consider «2 for large negative £ . We have 

P i / 2 


U 2 


, , pv- _ pff* , e -^le ,a -f 3 l e - W a -f 2 l 

< (i - c 2 ) b x 2V2 /_ de - 


' Vie* -pi 

•i-M/C dg 


£% - C l)Bl 2 V2 r l( 

IT •'■+<•/« ^|1 — g 3 | 


(A17) 


It is clear that the righthand side of (A17) goes to zero as £ — ► —00 and thus u 2 is 
approaches zero as £ — * — 00 . 

It is now appropriate to break up us(£) into two integrals: 


where 


Us(£) = u 3X (£) + u 3 2(£) 

««|f] = ~f° d(' «-»*<"-<*> [tf 0 (p| r |) + *,,(?,,)] 

7T Jf+u 


(Ais) 

(•419) 
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Now 


u 32 [ t] = £ ^(0 c-5 p («' 2 -« 2 ) [ifo(P|r|) + KoiPn)} {A20) 


“82 1 < 



1 

v/f’ - c 2 


< 


p »/* ^ 

2 V2B x C 2 



1 

\/£ J - £'* 


P 1 / 2 /— ri+C/M 1 

< — 2 V2 B, <M-Q- / o dq q- -j=^ (A21) 

It is easy to see that the right hand side of (A21) goes to zero as { -» -oo . Thus in 
this limit u 32 — > 0 . 

We now consider u 3i . First we have the known asymptotic property of the modified 
Bessel function Kq it is clear that for large enough argument, say larger than 10, one can 
choose constant C 4 , a pure number, so that 


It is convenient to break 



2, C 3 e-PH 

r 1 P»l* i r i»/» 


r c, ^ 

• 1 P3/2 ri 3/2 

up u 3 i into two more integrals: 


(A22) 

(A23) 


“31 — “311 + “312 


where 


and 


where 


and 


“311 


V* * U+, Vn 


“312 = 


VM 


m, a = 




vW 


\/M _ i 


(A24) 

(A25) 

(A26) 

(A27) 

( A28 ) 
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In view of (A4) and (A5), it is easy to show that in the range of integration each of 
!?(£, £') and Ei(£, £') are bounded above in absolute value by B 2 , for some 

constant B 2 and so the latter term in (A26) is 

< B 2 P- 1 ' 2 [° d? (£ 2 - f)" 1 ' 2 (A29) 

•'€+#» 

On substitution of £' = q£ into the integral in (A29), it is easily seen that the contribution 
from (A29) tends to 0 as £ — > -00 . Thus in the limit of £ -* — 00 , we are left only 
with the contribution from the 1st term in (A26) which is independent of xr and hence 
must be that from the Ivantsov solution. But it is known that the u is equal to A for the 
Ivantsov solution. Thus the limit of the first term in (A26) in the limit of £ — >• —00 must 
be A . Thus, the proof is complete. 
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